#ETHNICITY PLOTS
javanese =
  individual_level %>%
  filter(dist < 10000) %>%
  filter(dl01f == "A")

sundanese =
  individual_level %>%
  filter(dist < 10000) %>%
  filter(dl01f == "B")

model1_javanese <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing, data=javanese)
model2_javanese <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing, data=javanese)
model3_javanese <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing, data=javanese)
model4_javanese <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing, data=javanese)

model1_sundanese <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing, data=sundanese)
model2_sundanese <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing, data=sundanese)
model3_sundanese <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing, data=sundanese)
model4_sundanese <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing, data=sundanese)

model1_javanese <- coeftest(model1_javanese, vcov=cluster.vcov(model1_javanese, cluster = javanese$mkid14))
model2_javanese <- coeftest(model2_javanese, vcov=cluster.vcov(model2_javanese, cluster = javanese$mkid14))
model3_javanese <- coeftest(model3_javanese, vcov=cluster.vcov(model3_javanese, cluster = javanese$mkid14))
model4_javanese <- coeftest(model4_javanese, vcov=cluster.vcov(model4_javanese, cluster = javanese$mkid14))

model1_sundanese <- coeftest(model1_sundanese, vcov=cluster.vcov(model1_sundanese, cluster = sundanese$mkid14))
model2_sundanese <- coeftest(model2_sundanese, vcov=cluster.vcov(model2_sundanese, cluster = sundanese$mkid14))
model3_sundanese <- coeftest(model3_sundanese, vcov=cluster.vcov(model3_sundanese, cluster = sundanese$mkid14))
model4_sundanese <- coeftest(model4_sundanese, vcov=cluster.vcov(model4_sundanese, cluster = sundanese$mkid14))

model1_javanese <- tidy(model1_javanese) %>% mutate(outcome = "outcome1", level = "Javanese")
model2_javanese <- tidy(model2_javanese) %>% mutate(outcome = "outcome2", level = "Javanese")
model3_javanese <- tidy(model3_javanese) %>% mutate(outcome = "outcome3", level = "Javanese")
model4_javanese <- tidy(model4_javanese) %>% mutate(outcome = "outcome4", level = "Javanese")

model1_sundanese <- tidy(model1_sundanese) %>% mutate(outcome = "outcome1", level = "Sundanese")
model2_sundanese <- tidy(model2_sundanese) %>% mutate(outcome = "outcome2", level = "Sundanese")
model3_sundanese <- tidy(model3_sundanese) %>% mutate(outcome = "outcome3", level = "Sundanese")
model4_sundanese <- tidy(model4_sundanese) %>% mutate(outcome = "outcome4", level = "Sundanese")

ethnicity_plot <-
  bind_rows(model1_javanese, model2_javanese, model3_javanese, model4_javanese,
            model1_sundanese, model2_sundanese, model3_sundanese, model4_sundanese) %>%
  mutate(outcome_label = case_when(outcome == "outcome1" ~ "Place of worship",
                                   outcome == "outcome2" ~ "Live in village",
                                   outcome == "outcome3" ~ "Live in neighborhood",
                                   outcome == "outcome4" ~ "Rent room",
                                   TRUE ~ NA_character_)) %>%
  filter(term == "opium_legal") %>%
  mutate(Ethnicity = fct_rev(level)) %>%
  ggplot(aes(x=estimate, y = outcome_label, group = Ethnicity, shape = Ethnicity, color = Ethnicity)) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbarh(aes(xmin=estimate-1.67*std.error, xmax = estimate+1.67*std.error), position = position_dodge(width = 0.4), height = 0) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  xlab("Effect of Opium Concession System") +
  theme_bw() +
  scale_color_grey() +
  theme(axis.ticks.y.left = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title.y = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.title = element_blank(),
        axis.line.y.left = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))

ggsave("./_4_outputs/figures/figure_7.tiff", plot = ethnicity_plot,  width = 5, height = 3, dpi = 300, compression = "lzw")
